// 版权所有2016 Go作者。版权所有。
// 此源代码的使用受BSD样式
// 许可证的约束，该许可证可以在许可证文件中找到。

package big

import "math/rand"

// ProbablyPrime报告x是否可能是素数，
// 应用带有n个伪随机选择碱基的米勒-拉宾检验
// 以及贝利PSW检验。
// 
// 如果x是素数，则ProbablyPrime返回true。
// 如果随机选择x而不是素数，则ProbablyPrime可能返回false。
// 对于随机选择的非素数，返回true的概率最多为¼ⁿ. 
// 
// 对于小于2的输入，ProbablyPrime是100%准确的⁶⁴. 
// 有关错误概率的进一步讨论，请参见Menezes等人，《应用密码学手册》，1997年，第145-149页，
// 和FIPS 186-4附录F。
// 
// 可能prime不适合判断对手可能为愚弄测试而精心设计的素数。
// 
// 从Go 1.8开始，允许使用ProbablyPrime（0），并且仅应用贝利PSW测试。
// 在Go 1.8之前，ProbablyPrime只应用了Miller-Rabin测试，ProbablyPrime（0）惊慌失措。
func (x *Int) ProbablyPrime(n int) bool {
	// 注意上面的文档评论：
	// 更准确地说，贝利PSW测试使用
	// 超强卢卡斯测试作为其卢卡斯测试，但由于没有人知道
	// 如何在贝利PSW测试中区分任何卢卡斯测试
	// （它们在经验上都同样有效），这些细节不需要被记录或隐含地保证。评论中没有提到“The”Baillie PSW test 
	// 由于这种普遍的模糊性。

	if n < 0 {
		panic("negative n for ProbablyPrime")
	}
	if x.neg || len(x.abs) == 0 {
		return false
	}

	// primeBitMask记录的是<64的素数。
	const primeBitMask uint64 = 1<<2 | 1<<3 | 1<<5 | 1<<7 |
		1<<11 | 1<<13 | 1<<17 | 1<<19 | 1<<23 | 1<<29 | 1<<31 |
		1<<37 | 1<<41 | 1<<43 | 1<<47 | 1<<53 | 1<<59 | 1<<61

	w := x.abs[0]
	if len(x.abs) == 1 && w < 64 {
		return primeBitMask&(1<<w) != 0
	}

	if w&1 == 0 {
		return false // x是偶数
	}

	const primesA = 3 * 5 * 7 * 11 * 13 * 17 * 19 * 23 * 37
	const primesB = 29 * 31 * 41 * 43 * 47 * 53

	var rA, rB uint32
	switch _W {
	case 32:
		rA = uint32(x.abs.modW(primesA))
		rB = uint32(x.abs.modW(primesB))
	case 64:
		r := x.abs.modW((primesA * primesB) & _M)
		rA = uint32(r % primesA)
		rB = uint32(r % primesB)
	default:
		panic("math/big: invalid word size")
	}

	if rA%3 == 0 || rA%5 == 0 || rA%7 == 0 || rA%11 == 0 || rA%13 == 0 || rA%17 == 0 || rA%19 == 0 || rA%23 == 0 || rA%37 == 0 ||
		rB%29 == 0 || rB%31 == 0 || rB%41 == 0 || rB%43 == 0 || rB%47 == 0 || rB%53 == 0 {
		return false
	}

	return x.abs.probablyPrimeMillerRabin(n+1, true) && x.abs.probablyPrimeLucas()
}

// 概率PRIMILLERRABIN报告n是否通过
// 米勒-拉宾素数测试，使用伪随机选择的碱基。
// 如果force2为真，则其中一轮必须使用基数2。参见《应用密码学手册》，第139页，算法4.24。n已知为非零。确定q，k，使nm1=q<k 
func (n nat) probablyPrimeMillerRabin(reps int, force2 bool) bool {
	nm1 := nat(nil).sub(n, natOne)
	k := nm1.trailingZeroBits()
	q := nat(nil).shr(nm1, k)

	nm3 := nat(nil).sub(nm1, natTwo)
	rand := rand.New(rand.NewSource(int64(n[0])))

	var x, y, quotient nat
	nm3Len := nm3.bitLen()

NextRandom:
	for i := 0; i < reps; i++ {
		if i == reps-1 && force2 {
			x = x.set(natTwo)
		} else {
			x = x.random(rand, nm3, nm3Len)
			x = x.add(x, natTwo)
		}
		y = y.expNN(x, q, n)
		if y.cmp(natOne) == 0 || y.cmp(nm1) == 0 {
			continue
		}
		for j := uint(1); j < k; j++ {
			y = y.sqr(y)
			quotient, y = quotient.div(y, y, n)
			if y.cmp(nm1) == 0 {
				continue NextRandom
			}
			if y.cmp(natOne) == 0 {
				return false
			}
		}
		return false
	}

	return true
}

// 可能Primelucas报告n是否通过“几乎超强”Lucas概率素数检验，
// 使用Baillie OEIS参数选择。这对应于Jacobsen表格上的“AESLPSP”（以下链接）.
// 将此测试与以2为基数的米勒-拉宾/费马测试相结合，得出贝利PSW测试。
// 
// 参考文献：
// 
// 贝利和瓦格斯塔夫，“卢卡斯伪素数”，计算数学35（152），
// 1980年10月，第1391-1417页，特别是第1401页。
// https:
// 2000年3月，第873-891页。
// https:
// 
// Baillie，“超强卢卡斯伪素数”，OEIS A217719，https:
// 
// Jacobsen，“伪素数统计、表和数据“，http:
// 
// “贝利PSW素性测试”，http:
// （请注意，Nelly对“特强”的定义“测试给出了错误的雅各比条件，正如雅各布森指出的那样
// 
// Crandall and Pomerance，《素数：计算透视》，第二版，
// Springer，2005.
func (n nat) probablyPrimeLucas() bool {
	// Discard 0，1.
	if len(n) == 0 || n.cmp(natOne) == 0 {
		return false
	}
	// Two是唯一的偶数素数。
	// 已被调用方检查，但这里允许单独测试。
	if n[0]&1 == 0 {
		return n.cmp(natTwo) == 0
	}

	// Baillie是选择D，P，Q，
	// 的“方法C”，如https:
	// 尝试增加P≥ 这样，D=P²-4（so Q=1）
	// 直到雅可比（D，n）=1。
	// 在几次试验后，对非平方n的搜索有望成功。
	// 在超过预期的失败之后，检查n是否为平方
	// （这将导致雅可比（D，n）=1表示所有D不除以n）一般认为这是不可能的。如果我们得到一份报告，我们会要确切的数字
	p := Word(3)
	d := nat{1}
	intD := &Int{abs: d}
	intN := &Int{abs: n}
	for ; ; p++ {
		if p > 10000 {
			panic("math/big: internal error: cannot find (D/n) = -1 for " + intN.String())
		}
		d[0] = p*p - 4
		j := Jacobi(intD, intN)
		if j == -1 {
			break
		}
		if j == 0 {
			// d=p²-4=（p-2）（p+2）。
			// 如果（d/n）=0，则d与n共享一个素因子。
			// 由于循环以递增的p进行，并且以p-2==1开始，所以
			// 共享的素因子必须是p+2。
			// 如果p+2==n，则n为素数；否则，p+2是n的适当因子。
			return len(n) == 1 && n[0] == p+2
		}
		if p == 40 {
			// 如果n是正方形，我们将永远找不到（d/n）=-1。
			// 如果n是非正方形，我们平均只需几次尝试就可以找到d。
			// 尝试40次后，花点时间检查n是否确实是正方形。
			t1 = t1.sqrt(n)
			t1 = t1.sqr(t1)
			if t1.cmp(n) == 0 {
				return false
			}
		}
	}

	// “超强Lucas伪素数”的Grantham定义，在第页的Thm 2.3之后。876 
	// （上面的D，P，Q变成了Δ，b，1）：
	// 
	// /让U_n=U_n（b，1），V_n=V_n（b，1），Δ=b²-4。
	// 基b的一个超强Lucas伪素数是一个复合n=2^RS+Jacobi（Δ，n），
	// 其中s是奇数，gcd（n，2*Δ）=1，因此≡ 0模n和V_s≡ ±2模n，
	// 或（ii）V{2^t s}≡ 对于某些0，0 mod n≤ t<r-1。
	// 
	// 我们知道gcd（n，Δ）=1，否则我们会在上面找到雅可比（d，n）==0。
	// 我们知道gcd（n，2）=1，因为n是奇数。
	// 
	// 排列s=（n-Jacobi（Δ，n））/2^r=（n+1）/2^r。
	s := nat(nil).add(n, natOne)
	r := int(s.trailingZeroBits())
	s = s.shr(s, uint(r))
	nm2 := nat(nil).sub(n, natTwo) // n-2 

	// 我们采用“几乎超强”测试，检查上述条件
	// 美国除外≡ 0 mod n，这允许我们避免计算任何U_k值。雅各布森指出，也许我们应该做一个完整的超强测试：
	// “使用克兰德尔和波默朗斯方程3.13也可以恢复U_n:
	// U_n=D^-1（2V{n+1}-PV n）允许我们运行完整的超强测试
	// 以单个模块反转为代价。此计算在GMP中简单快速，
	// 因此，我们可以以与
	// 几乎相同的性能获得完整的超强测试。”计算Lucas序列V_s（b，1），其中：

	// 
	// V（0）=2 
	// V（1）=P 
	// V（k）=P V（k-1）-Q V（k-2）。
	// 
	// （记住，由于上面的方法C，P=b，Q=1。）
	// 
	// /通常V（k）=α^k+β^k，其中α和β是x²-Px+Q的根。
	// Crandall和Pomerance（P.147）观察0≤ J≤ k、 
	// 
	// V（j+k）=V（j）V（k）-V（k-j）。特别是要快速地将下标加倍：
	// 
	// V（2k）=V（k）²-2 
	// 
	// 因此，我们可以从k=0开始，在日志中累积到k=s₂（s） 步骤。英语短文：英语短文：英语短文：英语短文：英语短文：英语短文：英语短文：英语短文：英语短文：英语短文：英语短文：英语短文：英语短文：英语短文：英语短文：英语短文：英语短文：英语短文：英语短文：英语短文：英语短文：英语短文：英语短文：英语短文：英语短文：英语短文：英语短文：英语短文：英语短文：英语短文：英语短文：英语短文：英语短文：英语短文：英语短文：英语。英语短文：英语短文：英语短文：英语短文：英语短文：英语短文：英语短文：英语短文：英语短文：英语短文：英语短文：英语短文：英语短文：英语短文：英语短文：英语短文：英语短文：英语短文：英语短文：英语短文：英语短文：英语短文：英语短文：英语短文：英语短文：英语短文：英语短文：英语短文：英语短文：英语短文：英语短文：英语短文：英语短文：英语短文：英语短文：英语。检查V（s）≡ ±2（模块n）。
	natP := nat(nil).setWord(p)
	vk := nat(nil).setWord(2)
	vk1 := nat(nil).setWord(p)
	for i := int(s.bitLen()); i >= 0; i-- {
		if s.bit(uint(i)) != 0 {
			t1 = t1.mul(vk, vk1)
			t1 = t1.add(t1, n)
			t1 = t1.sub(t1, natP)
			t2, vk = t2.div(vk, t1, n)
			t1 = t1.sqr(vk1)
			t1 = t1.add(t1, nm2)
			t2, vk1 = t2.div(vk1, t1, n)
		} else {
			t1 = t1.mul(vk, vk1)
			t1 = t1.add(t1, n)
			t1 = t1.sub(t1, natP)
			t2, vk1 = t2.div(vk1, t1, n)
			t1 = t1.sqr(vk)
			t1 = t1.add(t1, nm2)
			t2, vk = t2.div(vk, t1, n)
		}
	}

	if vk.cmp(natTwo) == 0 || vk.cmp(nm2) == 0 {
		// 检查U（s）≡ 0根据Jacobsen的建议，应用Crandall和Pomerance方程3.13:⁻1（2v（k+1）-pv（k））
		// 
		// 因为我们正在检查U（k）=0，所以检查2v（k+1）==pv（k）模n，
		// 或pv（k）-2v（k+1）==0模n就足够了。
		t1 := t1.mul(vk, natP)
		t2 := t2.shl(vk1, 1)
		if t1.cmp(t2) < 0 {
			t1, t2 = t2, t1
		}
		t1 = t1.sub(t1, t2)
		t3 := vk1 // 偷取vk1，在下面不再需要
		vk1 = nil
		_ = vk1
		t2, t3 = t2.div(t3, t1, n)
		if len(t3) == 0 {
			return true
		}
	}

	// 检查V（2^t s）≡ 对于某些0，0 mod n≤ t<r-1。优化：V（k）=2是V（k'）=V（k）²-2的一个不动点，所以如果V（k）=2，我们可以停止：我们永远找不到未来的V（k）==0。
	for t := 0; t < r-1; t++ {
			return true
		}
		if len(vk) == 1 && vk[0] == 2 { // vk==2 
			return false
		}
		// k'=2k 
		// V（k'）=V（2k）=V（k）²-2 
		t1 = t1.sqr(vk)
		t1 = t1.sub(t1, natTwo)
		t2, vk = t2.div(vk, t1, n)
	}
	return false
}
